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ABSTRACT 

A  comprehensive  thermodynamic  framework  is  proposed  for  modeling  the  behavior  of  poly¬ 
meric  composites  exhibiting  viscoelasticity,  damage  and  permanent  deformation.  In  the 
present  formulation,  the  total  strain  is  decomposed  into  viscoelastic  and  permanent  com¬ 
ponents.  The  effect  of  damage  on  the  constitutive  behavior  is  introduced  through  direct 
coupling  between  damage  and  the  internal  state  variables  representing  deformation  of  the 
undamaged  (virgin)  material.  This  approach  accounts  for  time- depen  dent  damage  as  well 
as  damage- induced  changes  in  material  symmetry.  The  viscoelastic  part  of  the  proposed 
constitutive  model  is  first  formulated  for  a  general  case  of  damage  coupled  with  nonlinear 
viscoelasticity,  and  subsequently  specialized  to  coupling  with  linear  viscoelasticity.  Perma¬ 
nent  deformations  are  discussed  in  view  of  recent  experimental  findings.  Finally,  the  proposed 
constitutive  model  is  applied  to  experimental  data  for  swirl-mat  polymeric  composites. 
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1.  INTRODUCTION 


There  is  a  growing  interest  in  the  use  of  polymeric  composites  for  structural  applications. 
Some  reasons  for  such  growing  interest  are  the  significantly  lighter  weight  of  these  materials 
and  the  ease  with  which  they  can  be  manufactured.  In  all  applications  durability  plays  a 
crucial  role.  For  instance,  since  the  expected  life  of  an  automotive  vehicle  is  roughly  fifteen 
years,  polymeric  composites  intended  for  use  in  automotive  applications  are  required  to 
endure  for  the  same  length  of  time.  It  is  well-known  that  polymeric  composites  containing 
relatively  small  fractions  of  brittle  fibers  (roughly  30%  by  volume)  creep  viscoelastically  (e.g., 
Jerina  et  al.,  1979).  In  addition,  experimental  investigations  (e.g.,  Suvorova,  1985;  Schapery, 
1989;  Maire,  1992;  Smith  and  Weitsman,  1996)  indicate  that  polymeric  composites  undergo 
both  damage  and  permanent  deformation.  The  above  features  of  material  behavior  are  all 
life-limiting  factors  in  that  they  can  cause  excessive  dimensional  distortions  and  stiffness 
and  residual  strength  degradations.  The  aim  of  this  paper  is  to  establish  a  framework  for 
modeling  the  foregoing  features  of  material  behavior.  Such  a  framework  is  the  basis  for 
performing  durability  analyses,  which  play  a  major  role  in  engineering  design. 

Damage  observed  in  polymeric  composites  may  occur  in  the  form  of  microcracks  along 
fiber  bundles  and  also  microdebondings  along  the  fiber-matrix  interface  (e.g.  Smith  and 
Weitsman,  1996).  Up  to  the  present  time,  modeling  such  microcrack  damage  processes  has 
been  concerned  primarily  with  brittle  materials  (e.g.  rocks  and  concrete),  exhibiting  es¬ 
sentially  elastic  behavior,  and  metals  exhibiting  plastic  or  creep  behavior  (e.g.  Lemaitre, 
1984;  Chaboche,  1988;  Krajcinovic,  1989;  Ju,  1989;  Hansen  and  Schreyer,  1994).  A  com¬ 
mon  feature  in  these  works  is  that  the  effect  of  damage  on  material  behavior  is  introduced 
through  the  so-called  effective  stress.  The  notion  of  effective  stress  was  first  proposed  for 
one-dimensional  creep  problems  by  Kachanov  (1958),  who  introduced  the  effective  stress  as 
the  applied  stress  magnified  by  a  scalar  factor  representing  the  ratio  between  the  original 
(intact)  load  carrying  area  and  the  reduced  area  due  to  damage.  However,  it  is  well-known 
that  microcrack  damage  also  affects  the  material  symmetry.  To  account  for  this,  the  effective 
stress  was  subsequently  introduced  through  a  tensorial  mapping  of  the  applied  stress. 

Despite  the  progress  made  in  modeling  microcrack  damage  processes,  less  attention  has 
been  paid  to  modeling  such  processes  in  viscoelastic  materials  (e.g.  plastics  and  polymeric 
composites).  Notably,  Schapery  (1981,  1994,  1996)  established  the  basic  formulation  for 
viscoelastic  response  that  is  accompanied  by  microstructural  changes,  such  as  profuse  micro- 
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cracking.  These  formulations  account  for  viscoelasticity  by  means  of  a  set  of  scalar  valued 
internal  state  variables  similar  to  that  employed  in  this  article.  The  microstructural  changes 
are  represented  in  Schapery’s  work  by  means  of  another  set  of  internal  state  variables  whose 
evolutionary  laws  are  motivated  by  considerations  of  viscoelastic  fracture  mechanics.  An  ad¬ 
ditional  set  of  internal  state  variables  is  also  employed  by  Schapery  to  account  for  permanent 
(non-recoverable)  deformation.  It  may  also  be  mentioned  that  Weitsman  (1988)  attempted 
to  incorporate  the  effect  of  damage  by  coupling  viscoelasticity  and  damage  for  a  special  class 
of  linear  viscoelastic  materials. 

In  this  work,  the  formulation  of  coupled  viscoelastic  and  damage  response  is  patterned 
on  the  methodology  of  continuum  damage  mechanics  (e.g.  Lemaitre  and  Chaboche,  1985; 
Krajcinovic,  1989;  Lemaitre,  1992).  The  formulation  employs  several  of  Schapery’s  concepts 
as  well  as  the  thermodynamics  framework  outlined  by  Lubliner  (1972),  but  is  cast  within  the 
context  of  continuum  damage  mechanics.  This  format  is  correlated  with  damage  induced 
changes  in  material  symmetry  and  results  in  expressions  incorporating  effective  stresses. 

An  additional  issue  considered  in  this  work  is  the  permanent  deformation  that  accompa¬ 
nies  the  viscoelastic  damage  behavior  of  polymeric  composites.  In  general,  this  permanent 
deformation  increases  with  the  load  as  well  as  creep  duration  (e.g.  Tuttle  et  al.,  1995  ; 
Smith  and  Weitsman,  1996).  The  physical  origin  for  permanent  deformation  in  polymeric 
composites  arises  from  several  possible  sources.  For  instance,  the  ability  of  the  polymeric 
matrix  itself  to  undergo  irreversible  structural  changes  may  lead  to  non-recoverable  strains. 
Permanent  strains  could  also  be  due  to  frictional  sliding  following  fiber-matrix  interfacial 
microdebondings,  as  well  as  internal  surface  roughness  and  other  irregularities  resisting  mi¬ 
crocrack  closure  upon  unloading. 

In  Section  2  of  this  article,  we  present  a  general  thermodynamics  framework  that  ac¬ 
counts  for  viscoelasticity,  damage  and  permanent  deformation.  We  proceed  by  modeling  the 
coupling  between  viscoelasticity  and  damage  in  Section  3.  In  Section  4  permanent  defor¬ 
mation  is  considered  briefly  but  in  a  general  manner.  In  Section  5,  the  model  is  applied  to 
a  type  of  so-called  swirl-mat  polymeric  composites.  Section  6  concludes  with  some  general 
remarks  pertinent  to  the  present  work  and  also  contains  some  specific  observations  relevant 
to  the  material  considered  in  Section  5. 
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2.  GENERAL  THERMODYNAMICS  FRAMEWORK 


Consider  a  polymeric  composite  material  and  let  /yr  (r  —  1,2, . . . ,  R)  denote  R  scalar  valued 
internal  state  variables  representing  the  internal  degrees  of  freedom  of  molecular  motion  in 
the  polymeric  matrix.  Also,  let  (n  ( n  =  1,2,..., A)  denote  N  scalar  valued  internal  state 
variables  associated  with  permanent  deformation  l.  The  internal  state  variable  representing 
damage  is  taken  in  the  form  a  symmetric  fourth  rank  tensor  which  suffices  to  account 
for  general  damage-induced  changes  in  material  symmetry.  The  formulation  can  be  readily 
modified  to  accommodate  damage  variables  of  other  tensorial  ranks.  Throughout  this  paper 
the  subscripts  r  and  q  are  reserved  for  viscoelastic  quantities;  m  and  n  for  quantities  associ¬ 
ated  with  permanent  deformation;  and  a,  b,  c,  d,  i,j ,  Ic,  l,  p  and  s  are  associated  with  tensorial 
quantities  and  cover  the  range  1,2,3.  Also,  the  summation  convention  is  implied  over  the 
range  of  repeated  indices  unless  stated  otherwise. 

Viscoelasticity,  damage  and  permanent  deformation  are  all  forms  of  irreversible  thermo¬ 
dynamic  processes.  For  a  closed  system  and  small  strain  formulation,  the  entropy  production 
inequality  can  be  written  in  the  form  (Coleman  and  Gurtin,  1967) 

>  0,  (1) 

where  <f>  is  the  Gibb’s  free  energy  (per  unit  volume),  e,j  -  components  of  a  suitably  defined 
volume  average  infinitesimal  strain  tensor,  <x,_,  -  components  of  the  Cauchy  stress  tensor, 
S  -  entropy  (per  unit  volume),  T  -  temperature,  qt  -  components  of  the  heat  flux  vector, 
T,  =  dT/dxi  -  components  of  the  temperature  gradient,  and  a:,-  -  space  coordinates.  Also, 
in  (1)  the  overdot  signifies  differentiation  with  respect  to  time. 

Consider  a  Gibb’s  free  energy  of  the  form 


0  —  4*  (  )  Or  i  Cn  i  ^abcd  i  T )  . 


The  function  <f>  is  assmued  to  be  continuous  and  sufficiently  differentiable  with  respect  to 
its  arguments.  Considerations  of  the  entropy  production  inequality  in  (1)  give  the  familiar 
relations 


£ 


ij  — 


d<f> 

ton  ’ 


(2) 


^ere  we  adopt  the  variables  selected  by  Lubliner  (1990),  though  some  of  those  variables  may  in  fact  be 
components  of  the  second  rank  “back  stress”  tensor. 
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(3) 


5  = 


d£ 

dT  ’ 


and 

T 

rr7r  +  ^nC„  +  jki  <^ijki  —  ~"yr"  ^  0  5  (4) 

where  I\,  and  are  the  thermodynamic  forces  conjugate  to  the  internal  state  variables 
7r,  and  u>ijki,  respectively,  and  are  given  by 


and 


ijkl  — 


8<f) 


dui. 


'ijkl 


Guided  by  Lubliner  (1972),  we  assume  strain  decomposition  in  the  form 


(5) 

(6) 

(7) 


&ij  (&kl i  7j*  >  ^abcd y  T)  ~f*  &tj  ( Cft  i  ^abcdi  7°)  , 


(8) 


where  £?■  and  sfj  are  components  of  the  viscoelastic  and  permanent  strain  tensors,  respec¬ 
tively.  Since  <j>  is  a  continuous  function  of  atJ ,  then  it  follows  from  (2)  that 

dzij  _  d£ki 


dcrkl  d<7ij  ' 


Substitution  of  (8)  into  (9)  yields 


K  _  del, 


da u  da,.  ' 

Equation(lO)  implies  that  there  exists  a  scalar  valued  function,  say  </>v,  such  that 

d<j>v 


£ij  3a i 


(9) 


(10) 


(11) 


We  denote  by  <j>y  the  Gibb’s  free  energy  associated  with  the  viscoelastic  deformation. 
From  (2),  (8)  and  (11)  it  follows  that  <f>  can  be  written  in  the  general  form 


<j)  —  <p  a ,j  £jj  -|-  F{lr,  Cni^abcdi  T )  , 


(12) 


where  T  is  an  arbitrary  scalar  valued  function  of  its  arguments.  The  physical  significance  of 
J-  becomes  clear  when  we  consider  the  Helmholtz  free  energy  ip  using  the  relation  (Callen, 
1960) 

7r>  Cni  ^Oabcdi  1"}  =  ^P{^-iji1ri  Cni^abcdiT)  (Tij  £{j  .  (13) 
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Following  Lubliner  (1972),  we  assume  that  can  be  decomposed  in  the  form 


0  =  fl’''{£Vij,lr,^abcd,T)  +  4>P((n,Uabcd,T)  ,  (14) 

where  0V  and  0P  are  the  Helmholtz  free  energies  associated  with  the  viscoelastic  and  per¬ 
manent  deformations,  respectively.  Substituting  (12)  and  (14)  into  (13) 

<j>V(<Jij,'rT,Uabcd,T)+f{'rr,(n,Uabcd,T)  =  0V(e-i,7r,Watcrf,r)-<Ttj£-j-f0P(Cn,Wa6cd,r).  (15) 

The  Legendre  transform  of  0V  (Callen,  1960)  gives 


0  i^ij  1  TV?  ^abcdi (T{j  £jj  Q i^ij  ?  "Tr )  ^abcdi  T)  ,  (f  6) 

where  Q  is  a  scalar  valued  function  of  its  arguments.  From  (15)  and  (16),  it  is  clear  that 
Q  —  <f>v  and  T  —  0P;  whereby  T  is  independent  of  jr.  Thus  (16)  yields 


and  (12)  reads 


Since 


<t>V{(7ij,lr,Uabcd,T)  =  4>V(eVlV'fr,Uabcd,T)  -  CTtJ  £  £  , 

0  =  0V  -  atJ  ep-  +  0P  . 
d 


dc , , 


(-°ij  4  +  V,p)  =  4j  > 


equation  (18)  may  be  written  as 


(17) 

(18) 


4*  —  4*  {p’ij’i  Kt’i  ^abcch  T}  +  4*^ {&ij  i  Cni  ^abcdiT'j  , 


(19) 


where  <j>p  is  given  by 


and 


0P  =  0P  —  (Tij  £?j  . 


ep  =_Sf_ 

*  ’ 


In  view  of  (5)  and  (19),  the  thermodynamic  forces  Tr  are  given  by 

d<f>v 


Tr  = 


Similarly,  (6)  and  (19)  yield 


dir  ' 

<90p 

dCn  ’ 


(20) 

(21) 


(22) 

(23) 
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Also,  in  view  of  (7)  and  (19)  the  thermodynamic  forces  can  be  decomposed  in  the  form 


0 


ijkl  ~  *hjkl 


+nr. 


ijkl  ? 


(24) 


where 


df_ 

du>ijki  ’ 


and 


^ijkl  — 


df_ 

dujijki  ' 


Finally,  from  the  dissipation  inequality  (4),  we  have  the  following  requirements 


(25) 

(26) 


IV7r>0,  (27) 

rr7 r  +  >0,  (28) 

^nCn>0,  (29) 

and 

C„  +  f Ifjkl  ^ijkl  ^  0  •  (80) 

Inequality  (27)  should  always  be  satisfied  whenever  viscoelastic  deformation  occurs.  When 
viscoelastic  deformation  is  accompanied  by  damage  inequality  (28)  should  be  satisfied  as  well. 
Similarly,  inequalities  (29)  and  (30)  should  be  satisfied  whenever  permanent  deformation 
occurs  and  whenever  such  deformation  is  accompanied  by  damage,  respectively. 

This  establishes  the  general  thermodynamic  framework.  The  decompositions  (8)  and 
(19)  result  in  forms  where  the  coupling  between  viscoelasticity  and  damage,  and  between 
permanent  deformation  and  damage  are  split  into  two  separate  entities.  This  subject  is 
discussed  in  more  details  in  the  following  two  sections. 


3.  VISCOELASTICITY  COUPLED  WITH  DAMAGE 

For  the  sake  of  generality  and  to  encompass  a  majority  of  practical  applications,  the  vis¬ 
coelastic  constitutive  equations  will  be  first  formulated  for  a  general  case  of  nonlinear  vis¬ 
coelasticity  coupled  with  damage.  The  resulting  equations  will  subsequently  be  reduced  to 
linear  viscoelasticity  coupled  with  damage. 
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3.1.  General  Formulation 


The  formulation  will  be  first  obtained  for  the  case  of  fixed  stress  ertJ,  damage  u>ijki  and 
temperature  T.  The  results  will  then  be  extended  to  fluctuating  <r,j,  u>ijki  and  T. 

For  fixed  <r,j,  u)tjki  and  T,  an  irreversible  thermodynamic  process  is  triggered  in  the 
material,  which  prompts  the  viscoelastic  internal  state  variables  yr  to  drift  spontaneously 
toward  their  equilibrium  values  7*  =  uw*,T).  These  equilibrium  values  are  assumed 

to  be  continuous  and  sufficiently  differentiable  functions  of  their  arguments.  Assuming  that 
7r  and  7*  are  sufficiently  small,  a  Taylor  series  expansion  for  <f>v  about  7®  takes  the  form 

f  =  fi  +  \  (ir  -  it)  (7,  -  7;)  +  H  O  T.  ,  (31) 

where 

4>l  =  4 K((?ij,Uabcd,T ) 

is  the  value  of  (f>y  at  equilibrium,  and 


is  a  symmetric  matrix  assumed  to  be  constant.  In  the  above  relations,  the  subscript  “  e  ” 
implies  that  a  quantity  is  calculated  at  7r  =  7*  Vr.  Note  that  at  equilibrium  <f>y  is  minimum 
(Callen,  1960),  and  hence 

df_ 

dir 

and 

<f>rq6lrhq  >  0- 

Consequently,  there  is  no  linear  term  in  (31)  and  <f>vTq  is  a  positive  definite  matrix.  In  view 
of  (22),  it  is  clear  that  neglecting  the  higher  order  terms  (H.O.T.)  in  (31)  corresponds  to 
assuming  a  linear  relation  between  the  thermodynamic  forces  Tr  and  the  internal  state 
variables  iT.  Such  a  linear  relation  is  believed  to  be  valid  at  least  in  the  neighborhood  of 
equilibrium  (Prigogine,  1967).  In  Appendix  A,  we  also  show  that  neglecting  the  H.O.T.  in 
(31)  is  equivalent  to  neglecting  the  same  H.O.T.  in  the  internal  entropy  production.  It  should 
be  mentioned  that  an  expansion  similar  to  that  in  (31)  was  previously  used  by  Lubliner  (1972) 
in  a  strain  formulation  for  viscoelasticity. 
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Employing  the  usual  assumption  of  viscous-like  resistance  (Biot,  1954,  1955;  Schapery, 
1964,  1969),  let 

rr  =  ar,7,,  (32) 

where,  according  to  Onsager’s  principle  (Callen,  1960;  Fung,  1965;  deGroot  and  Mazur, 
1984),  arq  =  arq((Tij,ujabcd,T )  is  a  symmetric  matrix.  Substituting  (32)  into  inequality  (27) 

®rg  7r  7<7  ^  6  . 

Hence,  the  matrix  arq  is  positive  semi-definite. 

By  analogy  with  Biot  (1954,  1955)  and  Schapery  (1964,  1969),  we  assume  that  all  arq 
have  common  dependence  on  at] ,  u>ijki  and  T ;  namely 

Qrq  —  a((Tij,u>abcd)T^  •  aTq  ,  (33) 

where  a  is  a  positive  scalar  valued  function  of  its  arguments,  and  a°q  is  a  constant  symmetric 
and  positive  semi-definite  matrix.  Equations  (22),  (31)  and  (33)  yield 

+  «,%  =  «,  7,',  (34) 

where  £  is  a  reduced  time  defined  by  d£  =  dt/a.  Assuming  that  £  =  0  at  t  =  0,  then  for 
fixed  a ijki,  ujijki  and  T  one  has  £  =  t/a. 

Since  a°q  is  a  constant  symmetric  positive  semi-definite  matrix  and  <j>vTq  is  a  constant 
symmetric  positive  definite  matrix,  equation  (34)  can  be  rewritten  in  a  decoupled  form 
(Meirovitch,  1967)  as 

^r"7  +  ^ri  =  7r  (n°  SUm  OVer  r )  >  (35) 

d£ 

where  7r  is  a  transformed  set  of  internal  state  variables,  each  being  a  linear  combination 
of  the  original  internal  state  variables  7,.  The  parameters  7*  are  the  equilibrium  values 
corresponding  to  jr  and  are  obtained  from  7*  by  the  same  linear  transformation  as  that  for 
TV.  Also,  A®  and  are  constants  such  that  A°  >  0  and  >  0- 

For  fixed  stress,  damage  and  temperature,  the  solution  of  equation  (35)  is 

TV  =  7*  (l  —  e_^Tr)  (no  sum  over  r) ,  (36) 

where  rr  are  retardation  times  given  by 

A° 

tt  =  — ^  (no  sum  over  r )  .  (37) 

*  r 
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In  terms  of  the  transformed  internal  state  variables,  expansion  (31)  is  rewritten  as 


*'  =  «  +  \  £  (ir  -  7r')2  +  H.O.T. 


(38) 


The  viscoelastic  strain  can  now  be  obtained  by  substituting  (38)  into  (11)  bearing  in  mind 
that  7r,  and  hence  7,.,  are  to  be  kept  fixed  during  the  partial  differentiation  indicated  in  (11). 
Employing  (36)  we  then  obtain 


£ 


V 

ij 


d4>o 

dcij 


dAr 

dcr.j  ’ 


(39) 


where 

Ar  =  -  $*  (7r)2  (no  sum  over  r)  >  (40) 

and 

&  =  r.  +  £  a.  .  (41) 

r 

d>0  is  the  portion  of  the  Gibb’s  free  energy  associated  with  instantaneous  deformation. 

Recall  that  equation  (39)  was  obtained  for  fixed  stress,  damage  and  temperature.  For 
fluctuating  stress,  damage  and  temperature  the  strain  can  be  obtained  by  applying  the 
superposition  principle  to  expression  (39).  Following  Pipkin  (1986),  £  is  identified  as  an 
intrinsic  material  time  scale  for  the  viscoelastic  processes.  If  <7;j,  u>ijki  and  T  can  be  expressed 
in  terms  of  the  material  time  £  rather  than  physical  time  t,  then  dAr/d<Tij  becomes  a  function 
of  Straightforward  application  of  the  superposition  principle  (Pipkin,  1986)  yields 


£ 


V 

ij 


d4>o 

d<ji3 


+  j £  [l  -  e-*-™- 


if!, 


(42) 


where 


and  for  0  <  r  <  t 


dt' 

^  Jo  a[<Tij(P),Wa6c<i(P),?1(P)]  ’ 


(43) 

(44) 


f  _ E _ . 

Jo  a[<7,■j(t'),^a6ai(<')471(^,)] 

Noting  that  d£  =  dt/a  and  d£'  =  dr / a  the  variable  of  integration  in  (42)  can  be  reverted 
form  £'  to  r.  Hence 


Allowing  for  spatial  variations  of  stress,  damage  and  temperature,  the  total  derivative  d/dr 
in  (45)  is  replaced  by  a  partial  derivative,  holding  the  spatial  coordinates  fixed,  so  that 


e 


V 

ij 


d<t>o 

da. 


+ 


£?['- 


-K-O/Tr 


(46) 


(47) 


and  from  (43)  and  (44) 

<_<.=  /* _ dl _ . 

Jt  Cl  (  a%j  ,  ^abcd ,  T) 

The  first  term  on  the  right-hand  side  of  (46)  represents  the  instantaneous  part  of  the  strain, 
while  the  second  term  given  by  the  integral  represents  the  time-dependent  portion.  In 
Appendix  B,  it  is  shown  that  with  some  further  simplifying  assumptions  the  second  term 
can  be  considered  as  a  generalization  of  the  well-known  modified  superposition  integral  (e.g. 
Schapery,  1969). 

Using  (22),  (25),  (31),  (32)  and  the  transformed  internal  state  variables  7r,the  dissipation 
inequality  (28)  can  now  be  expressed  as 


d<t>0  /7r\ 

du>ijkl  “  V7r/ 


dAr 

du>ijki 


Vijki  >  0 . 


(48) 


The  complete  nonlinear  formulation  requires  an  expression  for  the  evolution  of  the  damage 
tensor  lOjjki  such  that  (48)  is  satisfied.  Such  expressions  can  be  formally  derived  from  ther¬ 
modynamic  considerations  through  the  introduction  of  a  damage  potential  (e.g.  Lemaitre 
and  Chaboche,  1985;  Lemaitre,  1992),  but  their  usefulness  seems  to  be  restricted  to  elastic 
response  with  damage.  Alternatively,  damage  evolution  relations  may  be  associated  with 
a  damage  function  akin  to  a  yield  function  in  plasticity  (Hansen  and  Schreyer,  1994).  An¬ 
other  approach,  which  is  guided  by  insights  from  viscoelastic  fracture  mechanics,  was  given 
by  Schapery  (1981,  1994,  1996).  This  approach  suggests  that  an  adequate  accounting  for 
damage  evolution  is  obtained  by  incorporating  damage  effects  only  within  the  instantaneous 
portion  of  ejj .  This  assumption  resulted  in  good  predictions  for  epoxy  resin  composites.  In 
practice,  the  evolution  equations  depend  on  the  particular  material  considered,  and  several 
such  expressions  can  be  found  in  the  literature  (e.g.  Krajcinovic,  1989).  A  form  of  these 
evolution  equations  will  be  considered  in  the  following  subsection. 
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3.2.  Linear  Viscoelasticity  Coupled  with  Damage 

Under  sufficiently  small  stress  levels,  the  viscoelastic  behavior  of  polymeric  composites  is 
linear.  In  this  section  the  foregoing  general  formulation  is  reduced  to  the  case  of  a  lin¬ 
ear  viscoelastic  behavior  coupled  with  damage.  First,  consideration  is  given  to  the  general 
case  of  anisotropy  with  respect  to  both  the  undamaged  material  behavior  and  the  ensuing 
damage  response.  Subsequently,  attention  is  given  to  the  simple  case  of  isotropic  material 
behavior  that  remains  isotropic  even  after  development  of  damage.  In  both  cases,  attention 
is  restricted  to  a  linear  viscoelastic  behavior  where  strain  depends  on  current  values  of  tem¬ 
perature  but  not  on  temperature  history.  Most  amorphous  polymers  exhibit  such  behavior 
(Ferry,  1980).  This  latter  restriction  simplifies  the  formulation  since  it  allows  us  to  discard 
coupling  between  damage  and  temperature  history.  The  extension  to  the  general  case  of 
temperature-history  dependent  strain  can  be  made  following  the  same  approach  adopted 
here. 


3.2.1.  Generally  anisotropic  behavior 


Consider  first  the  time-dependent  part  of  ejj.  Since  this  part  of  strain  is  assumed  to  be 
independent  of  temperature  history,  it  follows  that  all  7®,  and  hence  all  7®  and  Ar,  are 
independent  of  temperature.  Consequently 

Ar  —  Ar  (&ij ,  • 


Expanding  Ar  in  terms  of  stress  around  a  reference  state  (<77,  =  0,  T  =  Tr)  up  to  quadratic 
terms  to  retain  linearity  one  obtains 


d2Ar  \ 

dvijdcTki  J 


&ij  &ki  > 
0 


(49) 


where  the  subscript  0  implies  that  a  quantity  is  calculated  at  the  reference  state.  Note  that 
the  constant  term  in  (49)  vanishes  since  7®  =  0,  and  hence  7®  =  Ar  =  0,  at  the  reference 
state.  In  addition,  the  linear  term  in  (49)  is  discarded  since  it  corresponds  to  a  residual 
strain  at  the  reference  state,  which  is  disregarded  in  the  present  formulation. 

A  more  specific  functional  form  for  Ar  can  be  obtained  by  realizing  that  the  internal 
molecular  motions  represented  by  7r  occur  on  a  much  smaller  dimensional  scale  than  that  of 
damage  represented  by  u>ijki-  This  suggests  that  all  7*,  and  hence  all  7®  and  Ar,  are  likely  to 
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be  affected  by  damage  in  a  common  manner;  i.e.  they  have  common  dependence  on  ujabcd- 
Consequently,  we  can  write 


dcr a  da, 


ijV°kl )  o 


A -S' tips  Ppskl 


In  (50),  A Sijki  is  a  symmetric  fourth  rank  tensor  (A =  A Sfatj  =  A SjiU  =  A S![jlk),  and 
Pijki  =  Pijki{^abcd)  is  a  symmetric  fourth  rank  tensor  valued  function  of  the  damage  variable 
uJabcd  such  that 

at  UJabcd  =  0  *  Pijki  —  Pjkl  i  (hi) 


where 


lijkl  —  0  (&ik  fijl  T  &il  fijk)  i 


is  the  unit  fourth  rank  tensor  and  is  Kronecker  delta. 

In  view  of  the  above  hypothesis,  it  is  also  reasonable  to  assume  that  the  intrinsic  vis¬ 
coelastic  material  time  scale  £  is  unaffected  by  damage.  Moreover,  for  linear  viscoelastic 
behavior,  the  effect  of  stress  on  the  material  time  scale  is  negligible  (Ferry,  1980).  Thus,  the 
scalar  valued  function  a  in  (33)  takes  the  form 

a  =  ar(T) ,  (52) 

and  represents  the  usual  time-temperature  shift  factor  in  thermorheologically  simple  mate¬ 
rials  (Morland  and  Lee,  1960). 

Employing  (49),  (50)  and  (52),  the  viscoelastic  strain-stress  relation  (46)  can  be  written 


+  Jo_  A Sijki(t  -  {')  ^ ydr  , 


dai3 


where 


AS««(f)  =  £  (i  -  . 


is  the  time-dependent  compliance  tensor  of  the  undamaged  (virgin)  material, 

f  =  r  ?  =  r  _*i_ 

'  I ii  a\T(t')\  ’  5  Jo  a{T(t')\  ’ 


a ij  —  Pijki  &kl  • 
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Equation  (53)  suggests  that  the  viscoelastic  strain  e?  is  affected  by  damage  indirectly  through 
the  stress  d,j ,  which  is  the  applied  stress  <7,j  mapped  by  the  tensor  Pijki-  Following  common 
practice  in  continuum  damage  mechanics,  we  designate  <7,j  as  the  effective  stress. 

Consider  now  the  instantaneous  part  of  ejj.  The  equilibrium  viscoelastic  Gibb’s  free 
energy  <j>l  in  (38)  can  also  be  expanded  around  the  reference  state  in  the  form 


where 


«  =  -«&  Ar 


d(Tijd<Tkl 


j  &kl  ? 

0 


(57) 


\STdaij J0  ’ 

is  the  thermal  expansion  tensor,  AT  =  T  —  7r  is  temperature  excursion  from  its  reference 
value  Tr,  and  the  subscript  0  has  the  same  meaning  as  before.  Following  common  practice  in 
damage  mechanics  (e.g.  Lemaitre  and  Chaboche  1985;  Lemaitre,  1992),  we  further  assume 
that  ( d2<t>l/daijdcrki)o  depends  on  damage  in  the  separable  form 


\dat]  d<7kij  0 


S ijps  Qpskl  i 


(58) 


where  Qijki  =  Qijki(uabcd )  such  that  Qi]kl  =  Iijki  at  ojabcd  =  0,  and  Sfjkl  is  the  undamaged 
equilibrium  compliance  tensor.  Using  (41),  the  instantaneous  strain  takes  the  form 

-  A T  +  [Stjps  Qpskl  -  £  A S'jps  PPsk,}  <Tki ,  (59) 

0<7{j  T 


It  is  clear  that,  in  general,  the  instantaneous  part  of  strain  can  depend  on  damage  in 
a  manner  that  differs  from  that  of  the  time-dependent  part.  However,  for  simplicity  let  us 
assume  that  both  parts  have  the  same  dependence  on  damage  so  that 


Qijki  —  Pijki  • 


(60) 


From  (59),  the  instantaneous  strain  then  becomes 


d(Jij 


a°j  AT  + 


ijki  °ki 


where 

S?ir,  =  5'«  -  E  , 

r 

is  the  instantaneous  (elastic)  compliance  tensor  of  the  undamaged  material. 


(61) 

(62) 
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(63) 


The  total  strain  given  by  (53)  can  now  be  rewritten  in  the  compact  form 

r<  tn  daki 


4  =  a?-  AT  +  jf_  -  *')  -J^dr  , 


where 

5(ju(0  =  SJh  +  A Sijkl(( ) .  (64) 

is  the  overall  (instantaneous  and  time-dependent)  compliance  tensor  of  the  undamaged  ma¬ 
terial.  It  is  well-known  (e.g.  Fung,  1965)  that  both  and  Sfjkl  are  positive  definite.  From 
(62)  it  can  then  be  hypothesized  that  each  AS-*/  is  positive  semi-definite,  and  from  (54)  it 
follows  that  A Sijkiit)  is  also  positive  semi-definite,  as  it  should  be. 

In  the  remainder  of  this  work  we  consider  isothermal  conditions,  namely  AT  =  0.  The 
constitutive  equation  (63),  then  reads 

ejj  =  Jq_  Sijki(t  -  t)  .  (65) 

Also,  the  dissipation  inequality  (48)  becomes 

Ar^r  +  SijPs(t)  PpM  (T{j  Crkt  >  0  .  (66) 

r 

Since  the  first  term  on  the  left-hand  side  of  (66)  is  always  non-negative,  a  sufficient  condition 
to  satisfy  (66)  is 

Sijps{t)  Ppski  — >  positive  semi  —  definite ,  (67) 


where 


b  _  ^  Ppaki  . 

pskl  o  ^ abed  • 


dco, 


abed 


Thus,  the  functional  form  of  the  mapping  tensor  Pijki  is  restricted  by  the  two  requirements 
in  (51)  and  (67).  One  such  possible  functional  form  (e.g.  Ju,  1990)  is 


Pijki 


—  ( Pjkl  t^ijkl)  5 


(68) 


which  is  considered  to  be  a  generalization  to  the  scalar  form  used  by  Kachanov  (1958)  in 
one-dimensional  problems.  Clearly,  (68)  satisfies  (51),  and  the  choice  of  Wijki  should  be  such 
that  (67)  is  satisfied. 
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3.2.2.  Isotropic  behavior 


For  an  isotropic  undamaged  material  behavior,  the  instantaneous  compliance  in  (62)  takes 
the  form  (e.g.  Fung,  1965) 


=  +  5  (flo  -  !  j.)  Mu , 


(69) 


where  J0  and  Ba  are  the  instantaneous  shear  and  bulk  compliances,  respectively.  Similarly, 
the  time-dependent  compliance  in  (54)  takes  the  form 


A S«u(i)  =  j  A  J(t)  Iiju  +  i  [AB(()  -  |  A  J(<)]  4;  «u  , 


where,  analogously 


and 


A  J(t)  =  Y  A  Jr  (l  -  e-'"')  , 


(70) 

(71) 

(72) 


A B(t)  =  Y  A B,  (l  -  e-V*)  . 

r 

with  A  Jr  and  A Br  positive  constants.  The  overall  compliance  can  now  be  expressed  in  the 
form 

Sau(t)  =  i  J(t)  I*,,  +  i  [B(()  -  |  J(J)]  6,j 4, ,  (73) 

where 

J(t)  =  Jo  +  A J(t) , 

and 

B{t)  =  B0  +  A B(t) , 

are  the  shear  and  bulk  compliances,  respectively. 

Consider  the  circumstance  that  the  response  of  the  damaged  material  remains  isotropic. 
In  this  case,  damage  must  admit  an  isotropic  representation  that  does  not  change  the  original 
material  symmetry  properties.  One  possible  way  of  representing  isotropic  damage  is  to 
introduce  two  independent  scalar  variables  u>s  and  Ufj  that  characterize  the  influence  of 
damage  on  the  deviatoric  and  hydrostatic  strain  components,  respectively.  The  effective 


stress  ef ij  is  then  assumed  in  the  form 


CT,j  — 


U 


■  1  r 

' ~  o  /-.  \  °t j  i 


(1-075)  3(1  ~UH) 


(74) 


16 


where 


aij  ~  17  lJ 


-o'kk  l hj 


is  the  deviatoric  part  of  stress  and  &kk  is  the  trace  of  the  stress  tensor.  Also,  in  (74)  both 
and  u>h  are  in  general  functions  of  a  scalar  u>  representing  the  microcrack  density,  referred 
to  as  the  microcrack  concentration  parameter  (Kachanov, 1986). 

The  mapping  tensor  Pijki  can  be  evaluated  from  (56)  and  (74),  to  read 

1  „  (ufj  ~  US) 


Pijki  — 


lijkl  + 


(1  -us) 

Substitution  of  this  result  into  (68)  gives 


3(1—  U>5)  (1  —  uh) 


1 


Sij  Ski  ■ 


Uijki  =  us  Ujki  +  g  —  us)  8ij  hi  ■ 


(75) 


(76) 


Equation  (76)  takes  the  most  general  form  of  a  fourth  rank  isotropic  damage  tensor  ui,jki  (Ju, 
1990).  Upon  substitution  of  (73)  and  (74)  into  (65),  the  isothermal  strain-stress  relation  of 
an  isotropic  damaged  viscoelastic  response  takes  the  form 


dr 


(77) 


The  deviatoric  and  hydrostatic  parts  of  the  material  behavior  can  be  separated  from  (77) 
into  the  following  respective  forms 


2  </=/„'  <78) 

and 

where  £■)'  and  evkk  are  the  deviatoric  and  hydrostatic  strain  components,  respectively.  The 
effects  of  u>s  and  ujh  on  material  behavior  may  then  be  evaluated  from  uniaxial  creep  tests. 
For  instance,  let  <r,j  =  a  <5,i  8jl,  in  which  case  £52  =  e33>  then  (78)  and  (79)  yield 


and 


(80) 

(81) 
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Obviously,  equations  (80)  and  (81)  predict  a  stress  and  time-dependent  Poisson’s  ratio, 
v  =  i/(<7,  /  ).  In  principle,  these  equations  can  be  resolved  from  measurements  of  longitudinal 
and  transverse  strains. 

Isotropic  damage  thus  can  be  represented  by  the  two  scalar  functions  u>s  and  u>n.  These 
functions  are  subjected  to  two  restrictions.  Firstly,  recalling  equations  (51),  (68)  and  (76), 
they  must  vanish  in  the  absence  of  damage.  Secondly,  they  must  satisfy  the  requirement  in 
(67).  Using  (73)  for  Sijki(t)  and  (75)  for  Pijki,  the  left-hand  side  of  (67)  becomes 

Sijpa^  *>pM  =  (l  -I,s)2  “2^  {Iijkl  ~  3  6ij  Skl)  +  f  “4^  6ij  Skl '  ^ 

For  the  above  expression  to  be  positive  semi- definite,  it  suffices  that  both  and  Cofj  be 
non-negative.  This  implies  that  damage  can  only  increase  or  cease  to  increase  but  it  cannot 
decrease  (i.e.  no  microcrack  healing).  This  situation  encompasses  a  majority  of  practical 
applications  involving  damage. 

To  complete  the  isotropic  formulation,  evolution  equations  for  u>s  and  u>h  are  needed. 
Under  creep  loadings  with  increasing  damage,  such  equations  may  be  taken  in  the  well-known 
Kachanov- Rabotnov  form  (Kachanov,  1986) 


d>s 


As  (1  -  tc?s) 


(83) 


and 


uH  = 


1  Xh 


[Ah(1-uh)\ 


(84) 


In  equations  (83)  and  (84),  erf  and  a f*  are  equivalent  stresses  for  the  deviatoric  and  hydro¬ 
static  damage  modes,  respectively,  and  erf  and  considered  to  be  material  properties,  are 
the  corresponding  threshold  stresses  below  which  no  damage  develops.  Also,  xs,  Xh,  As, 
and  Ah  are  material  properties  characteristic  of  creep  damage,  which  are  typically  strongly 
dependent  on  temperature.  Finally,  (•)  denotes  the  McAuley  bracket,  i.e.  ( x )  =  x  if  x  >  0 
and  ( x )  =  0  if  x  <  0.  Guided  by  Leckie  and  Hayhurst  (1977),  the  equivalent  stresses  <rf  and 
<jH  may  be  taken  in  the  following  forms 


-f  = 


ij  -t  3 


(85) 


and 


<T?  =  +  (  1  -  /3)<Tk k  , 


(86) 
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where  <j\  is  the  maximum  principal  tensile  stress  and  /?(0</?<l)isa  material  property. 

For  monotonic  creep  loadings  the  stresses  a f  and  erf?  serve  as  thresholds  for  the  onset 
of  damage.  For  more  complex  loading  histories  involving,  for  instance,  complete  or  partial 
stress  removals,  the  onset  of  damage  is  usually  not  related  to  a  specific  threshold  stress.  In 
this  case,  the  maximum  strain  achieved  during  the  loading  history  may  serve  as  a  better 
indication  for  the  onset  of  new  damage  (Gurtin  and  Francis,  1981;  Suvorova,  1985).  For 
instance,  after  a  complete  or  a  partial  load  removal  damage  evolution  usually  ceases.  If 
afterwards  the  previously  achieved  maximum  strain  is  reached  by  increasing  the  applied 
load  or  by  creep  under  a  lower  load,  then  damage  would  start  to  evolve  again. 

Finally,  it  should  be  mentioned  that  a  simpler  representation,  where  isotropic  response  is 
retained  even  after  the  onset  of  damage,  can  be  obtained  as  a  special  case  of  (74)  by  setting 
ton  —  u>s  —  u.  Consequently,  damage  is  represented  by  the  scalar  u  so  that  u,jki  =  u>Iijki 
and  the  effective  stress  becomes 

**  =  ,  (87) 

1  —  u> 

This  effective  stress  coincides  with  the  original  form  proposed  by  Kachanov  (1958).  Under 
creep  loadings,  the  evolution  equation  for  u  can  again  be  taken  in  the  foregoing  Kachanov- 
Rabotnov  form  (Leckie  and  Hayhurst,  1977).  Note  that  in  contrast  with  the  elastic  case, 
where  scalar  damage  corresponds  to  a  constant  Poisson’s  ratio  (Ju,  1990),  in  the  viscoelastic 
case  Poisson’s  ratio  remains  both  time  and  stress  dependent. 


4.  PERMANENT  DEFORMATION 

The  viscoelastic  expressions  derived  in  the  foregoing  section  predict  strains  that  recover 
completely  after  load  removal  and  cannot  account  for  permanent  deformations  that  may 
persist  indefinitely.  As  mentioned  in  Section  1,  experimental  evidence  suggests  that  these 
permanent  strains  increase  with  load  during  monotonic  loading  as  well  as  with  creep  du¬ 
ration.  In  other  words,  permanent  strains  in  polymeric  composites  may  depend  on  both 
stress  magnitude  and  history.  The  foregoing  permanent  deformation  can  be  attributed  to 
several  deformation  mechanisms,  which  may  require  disparate  modeling  approaches.  In  the 
remainder  of  this  section  we  do  not  prescribe  a  specific  analytical  model  for  the  permanent 
deformation  but  provide  a  brief  review  of  some  of  the  possible  modeling  approaches.  The 
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thermodynamics  framework  formulated  in  Section  2  is  sufficiently  general  to  accommodate 
any  one  or  a  combination  of  these  modeling  methods. 

Since  permanent  strains  in  polymeric  composites  may  depend,  in  general,  on  both  stress 
and  stress  history,  £?■  in  (8)  should  be  further  decomposed  into  a  time-independent  compo¬ 
nent  and  a  time-dependent  component  £?• ,  namely 


£p  =  £?° 


+  £ 


(88) 


In  principle,  these  components  can  be  modeled,  respectively,  following  concepts  in  phe¬ 
nomenological  plasticity  and  viscoplasticity  (e.g.  Sun  and  Chen,  1989;  Ha  et  al.,  1991).  For 
instance,  experimental  investigations  (Smith  and  Weitsman,  1996)  reveal  that  under  uniaxial 
tensile  creep  loadings,  the  uniaxial  time-dependent  component  £pt  may  be  expressed  in  the 
form 

£pt  =  B  (<r  —  <7*)*  tp  ,  (89) 

where  <r*  is  a  threshold  stress  below  which  no  permanent  deformation  takes  place,  t  is  time, 
and  B,n  and  k  are  constants.  It  should  be  mentioned  that  an  expression  similar  to  that  in 
(89)  was  also  employed  by  Tuttle  et  al.  (1995)  for  the  time-dependent  permanent  deformation 
of  a  polymeric  composite.  Differentiating  (89)  with  respect  to  time  and  eliminating  t  in  favor 
of  £pt  and  (<r  —  a*)  we  obtain 


This  stress-strain  relation  is  in  the  well-known  power-law  form  of  viscoplasticity  (Lubliner, 
1990),  thereby  suggesting  that  ept  may  be  modeled  by  viscoplasicity  theory.  A  similar  ap¬ 
proach  can  be  taken  to  model  epo  using  time-independent  plasticity  instead  of  viscoplasticity. 

An  approach  to  model  plastic  deformation  in  viscoelastic  solids  has  also  been  proposed  by 
Naghdi  and  Murch  (1963).  This  modeling  effort  can  predict  the  onset  of  plastic  deformation 
under  general  circumstances  and  implicitly  accounts  for  the  decomposition  of  £p  in  (88). 
It  is  important  to  realize  that  adopting  phenomenological  plasticity  and  viscoplasticity  in 
modeling  the  permanent  deformation  associates  the  internal  state  variables  (n  in  Section 
2  with  internal  stresses  that  represent  isotropic  and  kinematic  hardening.  In  this  case, 
coupling  between  permanent  deformation  and  damage  can  be  modeled  following,  for  instance, 
approaches  by  Ju  (1989)  and  Hansen  and  Schreyer  (1994). 
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Physically  based  approaches  to  model  permanent  deformation  in  polymeric  materials 
have  also  been  proposed  (e.g.  Lee  et  al.,  1993;  Hasan  and  Boyce,  1995).  These  two  studies 
in  particular  focus  on  semi-crystalline  and  glassy  polymers,  respectively,  and  hence  could  be 
suitable  for  polymeric  composites  with  thermoplastic  matrices.  Also,  in  these  modeling  ap¬ 
proaches  it  is  assumed  that  the  permanent  deformation  is  associated  with  the  ability  of  the 
polymeric  matrix  to  undergo  irreversible  structural  changes  which,  correspondingly,  leads 
to  non- recoverable  strains.  In  polymeric  composites,  however,  permanent  deformation  can 
also  arise  from  other  mechanisms  such  as  frictional  sliding  following  fiber-matrix  interfacial 
microdebondings.  A  mechanistic  approach  for  modeling  such  phenomena  exists  for  unidi¬ 
rectional  fibers  embedded  in  an  elastic  matrix  (Hutchinson  and  Jensen,  1990).  No  similar 
models  exist  for  a  more  complex  fiber  architecture  (e.g.  randomly  oriented  fibers)  or  more 
complex  matrix  behavior  (e.g.  viscoelastic  matrix). 

Finally,  permanent  deformation  can  also  arise  following  microcracking.  This  is  attributed 
to  the  formation  of  microcrack  surface  roughness  and  other  irregularities,  due  to  fiber  ro¬ 
tations  and  breakage,  that  resist  microcrack  closure  upon  unloading  (e.g.  Schapery,  1989; 
Smith  and  Weitsman,  1996).  In  this  case,  the  internal  state  variables  £n  in  Section  2  can 
be  associated  with  the  internal  surface  roughness  and  the  internal  degrees  of  freedom  due  to 
fiber  rotations.  An  approach  similar  to  phenomenological  plasticity  can  also  be  adopted  to 
model  the  permanent  deformation  that  result  due  to  microcrack  damage  (Dragon  and  Mroz, 
1979). 


5.  APPLICATION  TO  SWIRL-MAT  POLYMERIC  COMPOSITES 

In  this  section,  the  foregoing  constitutive  model  will  be  utilized  to  predict  the  response  of 
a  swirl-mat  polymeric  composite.  The  objective  is  to  predict  the  creep-damage  behavior  of 
the  material  and  provide  guidelines  for  its  life-time  assessment  procedures. 

The  material  considered  consists  of  an  E-glass  fiber  preform  embedded  in  a  urethane 
matrix.  A  detailed  description  of  the  material  along  with  its  mechanical  response  is  given  by 
Smith  and  Weitsman  (1996).  Here,  we  only  observe  that  under  sufficiently  low  stress  levels 
the  material  exhibits  linear  viscoelastic  behavior,  while  under  applied  uniaxial  tensile  stresses 
that  exceed  a  threshold  level  of  ac  ~  36  MPa,  damage  in  the  form  of  profuse  microcracks 
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is  observed.  Also,  the  relative  magnitudes  of  the  permanent  strains  attained  during  creep 
tests  do  not  exceed  5%  of  the  total  strain.  Consequently,  from  the  practical  point  of  view, 
it  is  justified  to  neglect  the  permanent  deformation  for  the  present  purposes;  i.e.  e,j  ~  ejj. 

Experimental  data  for  the  material  under  consideration  are  available  for  uniaxial  tensile 
loadings.  Assuming  scalar  damage,  the  uniaxial  form  of  the  model  in  (65)  and  (87)  reads 

e  =  /„- D{t  ~ T)  i  (rb)  *  ’  (91) 

where  D(t)  is  the  uniaxial  compliance  of  the  undamaged  material.  Similarly  to  (64),  D(t) 
can  be  decomposed  into  instantaneous  D0  =  5°in  and  time-dependent  A D(t)  =  AS'im(t) 
parts,  namely 

D(t)  =  D0  +  AD(t) ,  (92) 


where,  in  view  of  (54), 

A D(t)  =  £  AS;m  (1  -  .  (93) 

r 

Upon  replacing  the  sum  of  the  discrete  spectrum  of  retardation  times  by  an  integral  of  a 
continuous  spectrum,  it  can  be  shown  (Pipkin,  1986)  that  A D{t)  may  be  expressed  as 


A  D(t)  =  D1tn, 


(94) 


where  D\  and  n(0  <  n  <  1)  are  constants,  which  can  be  readily  determined  from  uniax¬ 
ial  creep  tests  under  sufficiently  low  stress  levels  (i.e.  without  damage).  For  the  material 
at  hand,  Smith  and  Weitsman  (1996)  determined  D0  =  9.6  x  10-5MPa-1,  n  =  0.08,  and 
D\  =  2.8  x  10-5  MPa-1  /hour0  08.  A  variability  of  about  20%  is  observed  in  the  exper¬ 
imentally  determined  values  for  D0.  This  variability  is  attributed  to  randomness  in  the 
fiber  architecture  in  the  swirl-mat  as  well  as  to  manufacturing  induced  inhomogeneity.  The 
abovementioned  value  for  D0  is  therefore  taken  as  the  average  of  the  experimentally  recorded 
values. 

Similarly  to  (83)  and  (84),  the  evolution  equation  for  u>  is  taken  in  the  form 


u>  = 


(c  -  °c) 


X 


A[\  -u) 


(95) 


where  A  and  \  are  characteristic  creep  damage  properties.  It  should  be  noted  that  in  creep 
tests,  the  stress  is  typically  ramped  over  a  short  period  of  time  up  to  the  prescribed  level. 
During  the  ramp  loading  stage  the  behavior  is  essentially  linear  elastic  until  the  stress  exceeds 
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the  threshold  value  ac  when  damage  starts  to  occur.  Thus,  for  tests  with  a  >  crc,  creep  occurs 
under  an  initial  value  of  damage  u0  that  accumulates  during  the  load  ramping  stage  and 
depends  on  the  stress  level  a.  For  creep  tests  with  a  >  <xc,  solution  of  (95)  using  the  initial 
condition  u>  =  u)0  at  t  =  0  reads 

1  1  (\  iy1/{1+x) 

1  —  U>  1  —  U>0  V  tcJ 

where  tc  is  a  normalizing  constant  given  by 

_  (i^  . 

1+X  \(T-(JCJ 

Denoting  the  value  of  u  corresponding  to  a  complete  failure  of  the  material  by  u>f,  then  the 
time  to  failure  t[  can  be  obtained  from  (96)  in  the  form 


t( 


(98) 


From  (91)  and  (96)  it  follows  that  a  theoretical  value  for  u>  at  complete  failure  is  uc  =  1 
and  occurs  at  t  =  tc.  In  practice,  however,  materials  typically  fail  at  values  of  LO{  <  1  (e.g. 
Lemaitre  and  Chaboche,  1985,  p.  364)  corresponding  to  t{  <  tc. 

To  complete  the  uniaxial  characterization  of  the  material,  values  for  cj0,  u>{,  x  and  A  are 
needed.  If  w0  is  known  then  creep-rupture  tests  can  be  performed  at  different  stress  levels  to 
experimentally  determine  t{  as  a  function  of  stress,  and  subsequently  determine  the  values 
of  W{,  x  and  A. 

Focusing  attention  on  the  instantaneous  damage  u0  that  occurs  during  the  rapid  ramp 
loading  to  the  stress  level  <r,  we  note  that  in  the  case  of  elastic  materials  undergoing  damage, 
the  ratio  of  the  unloading  compliance  to  the  loading  compliance  provides  a  measure  for  the 
level  of  damage.  In  particular,  it  is  well-known  (Kachanov,  1986)  that 

1  iff 

1  -  u;0  D0  ’ 

where  £)“  is  the  unloading  compliance.  To  determine  ui0  “spike”  tests,  each  consisting  of 
a  constant  stress  rate  loading- unloading  cycle,  were  performed  up  to  the  stress  levels  a  = 
55,  69,  and  83  MPa  corresponding,  respectively,  to  approximately  40,  50,  and  60%  of  the 
ultimate  tensile  strength.  The  rate  of  loading  used  was  a  «  5.6  MPa/s,  which  is  sufficiently 
high  to  keep  the  viscoelastic  effects  at  a  minimal  level.  From  these  tests,  values  for  D™  were 
determined.  Values  for  D0  were  also  determined  to  eliminate  the  effect  of  sample  to  sample 
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Figure  1:  Ratio  between  the  uniaxial  unloading  and  loading  compliances  as  a  function  of 
stress 


variability  in  D0.  Figure  1,  which  shows  the  results  of  these  three  tests,  clearly  suggests  a 
linear  correlation  between  D"/D0  and  a.  Since  for  a  <  ac  no  damage  occurs,  then  using 
(99)  it  is  possible  to  write 


_1 _ 

1  -  u0 


=  1  + 


(°~  ~  gc) 

c 


(100) 


where  C  is  a  normalizing  constant  determined  from  Figure  1  as  C  —  510  MPa. 

With  the  above  expression  for  u0  at  hand,  it  is  now  possible  to  evaluate  o»f,  A,  and  \ 
from  time-to-failure  data.  Preliminary  creep-rupture  tests  at  temperature  T  =  75°F  and 
50%  relative  humidity  (Corum,  1995)  are  shown  in  Figure  2  .  Clearly,  these  data  exhibit 
a  large  amount  of  scatter  which,  again,  is  attributable  to  inhomogeneity  and  randomness 
of  the  swirl-mat  polymeric  composite  material.  Nevertheless,  the  data  in  Figure  2  can  be 
used  to  provide  some  life-time  estimates  for  the  material.  The  constants  u>t,  A,  and  x  were 
determined  by  fitting  the  experimental  data  to  the  expression  for  t[  in  (98)  using  (100)  for 
uj0.  Based  on  the  best  fit  curve,  depicted  in  Figure  2  by  the  solid  line,  the  following  values 
were  determined 


u>[  =  0.67  ,  x  =  7.1  ,  A  =  260  MPa.hour1/71  .  (101) 


In  view  of  the  large  scatter  in  the  creep-rupture  data,  the  above  values  should  be  considered 
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Figure  2:  Time  to  failure  as  a  function  of  stress  in  uniaxial  creep-rupture  tests 


as  preliminary  estimates  rather  than  conclusive  material  properties. 

In  addition  to  the  aforementioned  inherent  variability  in  properties  of  the  swirl-mat 
composite,  a  factor  that  significantly  contributes  to  the  scatter  in  creep-rupture  data  is 
the  high  sensitivity  of  the  behavior  of  the  considered  material  to  fluctuations  in  the  ambient 
environment  (temperature  and  relative  humidity).  Creep-rupture  tests  typically  require  long 
durations,  where  uncontrollable  fluctuations  in  the  environment  accumulate  and  can  affect 
the  results  significantly.  It  is  apparent  that  the  large  scatter  in  the  limited  experimental  data 
in  Figure  2  undermines  the  reliability  of  predictions  based  upon  the  values  listed  in  (101). 
We  therefore  restrict  attention  only  to  times  that  are  short  in  comparison  with  the  time  to 
failure  if,  so  that  the  effect  of  environmental  fluctuations  can  be  neglected. 

Short  term  creep  tests  (approximately  170  hours)  were  performed  under  stress  levels  of 
55,  69,  and  83  MPa  at  temperature  T  =  75°F  and  50%  relative  humidity.  The  above  stresses 
were  ramped  at  the  same  rate  as  that  for  the  tests  in  Figure  1,  so  that  the  expression  for 
uj0  in  (100)  remains  applicable.  For  creep  under  constant  stress  cr,  substitution  of  (92),  (94) 
and  (96)  into  (91)  yields 


e  =  Da 


(i  r^j +  D'  (T^j  L- (t  “ T)”  Tt  (* "  0 


-i/(i+x) 


dr 


(102) 
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Upon  integrating  by  parts  and  changing  the  integration  variable  so  that 


V  =  r/t, 


equation  (102)  becomes 


Dn 


cr  r1  ,  {  t  x-VU+x) 

+  nD'e  0^7)1  (1-,r  d’1' 


(103) 


(i-u-r  * 

The  solution  of  the  integral  in  (103)  can  be  obtained  in  a  closed  form  (Gradsteyn  and  Ryzhik, 
1965,  Sec.  3.211,  p.  287  and  Sec.  9.182-10,  p.  1055),  so  that 


e(t)  =  D0 


+  Dx  tn 


F(l,p,l  +  n;t/tc) , 


(104) 


(1-w)  (l-^o) 

where  F(l,p,l  +  n;  t/tc)  is  a  hypergeometric  series  that  converges  for  all  t  <  tc  and  is  given 

by 


F(l,  p,l  -f  n;  t/fc)  —  1  -f 


(n  +  1) 


T>0 


+ 


P&  +  1) 


(n  +  l)(n  +  2) 


& 


+ 


and  p  is  given  by 


P(P  +  1){P  +  2) 

(n  +  l)(n  +  2)(n  +  3) 

1 


(0‘ 


+ 


P  = 


1  +  X 


Experimental  data  along  with  model  predictions  are  given  in  Figure  3.  The  experimental 
data  are  depicted  by  the  scattered  symbols  while  model  predictions  are  represented  by  the 
solid  lines.  Note  the  good  agreement  between  model  predictions  and  experimental  results. 

Using  (98),  (100)  and  the  values  for  the  creep-rupture  constants  in  (101),  the  times 
to  failure  are  estimated  as  t{  1.07  x  107, 1.72  x  105  and  1.14  x  104  hours  for  the  stress 
levels  55,  69,  and  83  MPa,  respectively.  Obviously  a  creep  duration  of  170  hours  represents 
only  a  very  small  fraction  of  the  life  time,  during  which  the  damage  u>  remains  essentially 
constant  at  its  initial  value  uQ,  as  can  be  verified  from  (96).  This  is  in  agreement  with  the 
experimental  observation  (Smith  and  Weitsman,  1996)  that  after  creep  periods  ranging  from 
0.5  to  170  hours  the  unloading  compliance  remains  essentially  unchanged,  implying  that 
damage  remains  practically  constant.  Thus  for  short-time  creep-damage  behavior  a  good 
approximation  to  the  constitutive  equation  (91)  with  constant  a  is 

1 


1  —  UJr 


D(t)  a . 


(106) 
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Figure  3:  Comparison  between  model  and  experiment  for  short-term  creep 

Indeed,  using  (100),  the  difference  between  this  approximation  and  exact  calculation  of  the 
integral  in  (91)  is  quite  negligible  for  the  creep  durations  shown  in  Figure  3. 

For  illustration  purposes,  model  prediction  for  long-term  creep  behavior  up  to  the  failure 
time  tf,  as  obtained  from  (91),  for  a  =  83  MPa  is  shown  in  Figure  4.  This  figure  demonstrates 
that,  similarly  to  metal,  the  material  at  hand  exhibits  a  significant  amount  of  tertiary  creep 
prior  to  failure.  Features  of  such  behavior  were  observed  in  creep-rupture  tests,  especially 
at  high  temperatures  (Corum,  1995).  The  fact  that  tertiary  creep  can  be  significant  for  the 
material  considered  herein  is  also  consistent  with  the  relatively  large  value  found  for  u j. 
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6.  CONCLUDING  REMARKS 

In  this  article  a  thermodynamically  consistent  framework  is  proposed  for  modeling  different 
types  of  deformation  of  polymeric  composites.  The  framework  accommodates  viscoelasticity, 
damage  in  the  form  of  distributed  microcracks,  and  both  time  independent  and  dependent 
permanent  deformation.  Several  tensorial  ranks  for  the  damage  variable  can  be  used  in  the 
context  of  the  present  formulation.  In  this  article,  however,  damage  is  represented  by  a 
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Figure  4:  Model  predictions  for  creep-damage  behavior  up  to  failure 


general  fourth  rank  tensor.  Such  a  representation  is  capable  of  simulating  general  changes 
in  material  symmetry  properties  due  to  microcracking-induced  damage. 

Upon  establishing  the  thermodynamics  framework,  a  comprehensive  approach  was  pro¬ 
posed  for  modeling  the  coupling  between  viscoelasticity  and  damage.  The  formulation  was 
first  established  for  a  nonlinear  viscoelastic  behavior  coupled  with  damage.  It  is  shown  in 
Appendix  B  that,  in  the  absence  of  damage,  the  well-known  modified  superposition  integral 
is  a  special  case  of  the  present  formulation.  The  formulation  is  subsequently  reduced  to 
linear  viscoelastic  behavior  coupled  with  damage.  In  this  case,  the  effect  of  damage  can  be 
incorporated  into  the  constitutive  equations  in  a  form  consistent  with  the  well-known  effec¬ 
tive  stress  concept.  For  linear  viscoelastic  behavior  coupled  with  damage,  two  special  cases 
of  the  generally  anisotropic  damage  behavior  were  considered.  In  the  first  of  these  two  cases 
damage  is  represented  by  two  scalar  functions  tog  and  cjh  and  was  referred  to  as  isotropic 
damage.  In  the  second  case  damage  was  represented  by  a  single  scalar  and  was  referred  to 
as  scalar  damage. 

Using  scalar  damage,  a  uniaxial  form  of  the  model  for  linear  viscoelasticity  coupled  with 
damage  was  applied  to  a  swirl-mat  polymeric  composite.  When  the  deformation  model 
was  complemented  by  a  Kachanov-Rabotnov  form  of  damage  evolution,  an  experimentally 
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observed  tertiary  creep  stage  could  be  predicted  by  the  current  formulation.  This  affords 
a  convenient  means  for  life-time  assessment  of  structural  components.  It  is  noted  that  the 
estimated  strain  to  failure  of  about  2.3%  (Figure  4)  is  consistent  with  experimental  data. 
This  justifies  the  small  strain  formulation  adopted  in  this  work  which,  for  the  case  of  swirl- 
mat  polymeric  composites,  remains  valid  up  to  failure. 

It  is  interesting  to  note  that  a  special  sub-case  of  the  present  formulation  has  certain  simi¬ 
larities  with  a  somewhat  reduced  form  of  Schapery’s  non-linear  viscoelastic  model  (Schapery, 
1969).  To  show  this,  consider  the  general  strain-stress  relation  given  by  (46)  and  (47).  Adopt¬ 
ing  the  functional  forms  in  (49),  (50),  (57)  and  (58),  and  allowing  for  Qijki  #  Pijkh  then  the 
strain-stress  relation  (53)  remains  valid  with  the  instantaneous  strain  given  by  (59).  For 
simplicity,  consider  isothermal  conditions  AT  =  0.  Restricting  attention  to  scalar  damage 
and  assuming  that  the  material  time  scale  £  is  unaffected  by  damage,  then  similarly  to  (91) 
the  viscoelastic  uniaxial  strain-stress  relation  can  be  written  as 

ev  =  f0  D0 a  +  A£>(£  -  (')  a)  dr  ,  (106) 

where  f0  and  f<i  are  scalar  functions  of  the  damage  variable  and 

t=r 

Jo 

If  the  tensor  Ptjk,  is  taken  in  the  form  given  by  (68),  then  /2  =  1/(1  -  oj).  However,  since 
Qijki  ^  Pijki  in  (106),  f0  and  /2  can  be  disparate  functions  of  u.  Since  the  evolution  equation 
for  u  can  be  written  in  a  general  form  as  u  =  c 5(<r,u;),  then  u  can  in  principle  be  expressed 
as  a  function  of  a.  Consequently,  fQ  and  /2  can  be  considered  as  compounded  functions  of 
a  and  t. 

Schapery’s  model  (1969),  for  isothermal  conditions  and  applied  uniaxial  stress,  reads 

e  =  g0  D0  a  +  gx  J^_  A D(£  -  £’)  -£-( g2  a)  dr  ,  (107) 

where  g0,gi  and  gi  are  scalar  functions  of  stress  (but  not  of  time).  Comparison  between 
(106)  and  (107)  shows  that  the  two  forms  are  similar  with  the  exception  that  in  (106)  the 
counterpart  of  g\  can  be  thought  of  as  a  constant  lumped  within  /2.  It  is  worth  noting  that 
in  the  present  formulation  the  physical  origin  of  the  scalar  functions  /„  and  /2  stems  from 
damage.  By  contrast,  in  Schapery’s  formulation  the  functions  gQ,  g\  and  </2  arise  from  a  formal 
thermodynamic  formulation.  It  should  also  be  noted  that  unlike  Schapery’s  formulation,  in 


dt'  ,  _  r  dP 

a[<r( t')]  ’  Jo  a[cr(t')] 
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the  current  work  the  Gibb’s  free  energy  function  is  expanded  about  the  equilibrium  values 
7*  (equation  31),  which  leads  directly  to  expression  (50)  for  a  damage  modified  compliance 
familiar  in  the  context  of  continuum  damage  mechanics. 

While  emphasis  in  this  work  was  placed  on  the  viscoelastic  part  of  the  deformation, 
several  deformation  mechanisms  and  modeling  approaches  to  the  permanent  deformation 
were  also  discussed  and  can  be  readily  accommodated  within  the  proposed  thermodynamics 
framework.  For  simple  loading  cases,  e.g.  monotnoic  and  creep  loadings,  applied  to  the 
swirl-mat  polymeric  composite  considered  in  Section  5,  experimental  observations  indicate 
that  permanent  strains  are  significantly  lower  than  their  viscoelastic  counterparts  and,  hence, 
permanent  deformation  can  justifiably  be  neglected.  It  should,  however,  be  noted  that  under 
cyclic  loading  conditions,  and  for  a  sufficiently  large  number  of  cycles,  permanent  deformation 
may  play  an  important  role  in  the  overall  deformation. 
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In  this  appendix  it  is  shown  that  neglecting  the  higher  order  terms  (H.O.T.)  in  (31)  is  equiv¬ 
alent  to  neglecting  the  same  H.O.T.  in  the  internal  entropy  production.  To  this  end,  and  to 
make  the  point  in  the  simplest  possible  way,  consider  the  circumstance  where  no  permanent 
deformation  occurs.  Thus  <f)y  =  <j),  t/>v  =  ip  and  we  have  (e.g.  Lubliner,  1972) 

(f>v  =  u  —  aij  £{j  —  ST ,  (A.l) 

where  u  is  the  internal  energy  per  unit  volume.  To  further  simplify  matters,  consider  the 
case  of  fixed  stress,  damage  and  temperature,  then  at  equilibrium 

<f>l  =  Ue-  (Tij  (£,j)e  -  SeT  ,  (A.2) 

where,  again,  the  subscript  “e”  implies  that  a  quantity  is  evaluated  at  7r  =  7*.  Subtracting 
(A.2)  from  (A.l)  one  gets 

<t>v  -  <f>l  =  u  -  ue  -  ai}  [eij  -  (£tj)e]  -  T  AS  ,  (A. 3) 


where 


AS  =  S  -  Se, 


is  the  change  in  entropy  from  its  equilibrium  value,  which  can  be  decomposed  (Prigogine, 
1967)  in  the  form 
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AS  =  AeS  +  , 


(A.4) 


where  A eS  and  A;. S'  are  the  changes  (from  equilibrium)  in  external  entropy  supply  and 
internal  entropy  production,  respectively. 

The  first  law  of  thermodynamics  can  be  written  in  the  form  (Fung,  1965) 

du  —  Oij  deXJ  -f  T  deS ,  (A. 5) 

where  deS  is  an  infinitesimal  change  in  external  entropy  supply.  Integrating  the  above  dif¬ 
ferential  form  from  equilibrium  to  an  arbitrary  current  state  we  get 

u  —  ue  —  crij  [eij  -  (e,j)e]  +  T  A eS .  (A.6) 

Substituting  (A.6)  into  (A.4)  and  using  (A. 3)  we  obtain 

(A.7) 


For  fixed  stress,  damage,  and  temperature;  the  internal  entropy  production  takes  the 
form  (Prigogine,  1967) 

d,S  =  ^  ■  (A.8) 

Following  Prigogine  (1967),  expand  Tt  in  terms  of  7 r  around  the  equilibrium  values  7*,  not¬ 
ing  that  Tr  =  0  at  7r  =  7*  and  keeping  up  to  linear  terms  only,  to  get 

r.=  (|0  (7,-7,')  +  H.0.T.  (A.9) 

Substituting  (A.9)  into  (A.8)  and  integrating  from  equilibrium  to  the  current  state  yields 
(Prigogine,  1967,  pp.  47,  48) 

T  A;S  =  i  (jSj  hr  -  it)  (7,  -  7,')  +  H.O.T.  (A. 10) 

Clearly,  substitution  of  (A. 10)  into  (A.7)  and  using  (22)  gives  the  expansion  in  (31),  which 
establishes  the  stated  equivalence. 


APPENDIX  B 

In  this  Appendix,  we  show  that  a  special  form  of  (46)  can  be  viewed  as  a  generalization  of 
the  modified  superposition  integral.  To  this  end,  consider  the  special  case  in  which  all  7®, 
and  hence  all  7®  and  Ar,  have  common  dependence  on  stress,  damage  and  temperature,  i.e. 
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^rfaij  i  ^abcdt  T)  Cr  A(<7,j,Cc>0(,cj,  T)  Vr  ,  (B.l) 

where  <v  is  a  scaling  constant  and  A  is  a  scalar  function  of  its  arguments.  Denoting  the 
time-dependent  part  of  the  strain  in  (46)  by  and  using  (B.l),  we  can  write 


where 


(B.2) 


HO  =  5>  (l  -  e-i,r’)  ■  (B-3) 

r 

Considering  the  isothermal  undamaged  material  behavior  and  assuming  that  a  is  indepen¬ 
dent  of  ,  then  (B.2)  becomes 


e 


c 


h(t  —  t) 


d_ 

dr 


(B.4) 


and  A  becomes  a  function  of  the  stress  invariants  only. 

The  meaning  of  the  function  h(t)  can  be  established  by  considering  the  response  under 
sufficiently  low  stress  levels.  In  this  case  the  viscoelastic  behavior  is  linear  in  stress,  hence 
A  ^  cr2,  namely 


(B.5) 


where  A 5°  w  is  a  symmetric  fourth  rank  tensor.  Substitution  of  (B.5)  into  (B.4)  gives 

=  Jo_  A Sijki(t  -  t )  ,  (B.6) 

where 


A  Siiu(i)  =  h(t)AS°ju.  (B.7) 

Since,  by  hypothesis,  the  retardation  times  rT  are  independent  of  stress,  then  ASijki(t)  is  the 
time-dependent  compliance  tensor.  All  components  of  this  compliance  tensor  have  common 
time  dependence  expressed  by  h(t).  It  follows  that  /i(£),  which  is  analogous  to  AD(<)  in 
(94),  is  the  linear  uniaxial  creep  compliance  of  the  material  and  can  be  determined  from  a 
uniaxial  creep  test  under  sufficiently  low  stress  levels. 

With  the  above  interpretation  for  h(t)  and  assuming  that  a  is  independent  of  <7jj,  i.e. 
a  =  a^abcd^T),  it  is  clear  that  (B.2)  can  be  considered  as  a  multiaxial  generalization  of 
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the  modified  superposition  integral  which  also  accounts  for  damage  and  thermal  effects. 
This  interpretation  for  (B.2)  is,  however,  restricted  to  materials  characterized  by  common 
time-dependence  in  all  components  of  the  creep  compliance  tensor  in  the  linear  range  of 
behavior.  For  isotropic  behavior,  this  corresponds  to  the  case  of  materials  exhibiting  a 
constant  Poisson’s  ratio,  which  applies  for  many  polymers  (Ferry,  1980). 
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